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We study the effects of perturbative reheating on the evolution of the curvature perturbation 
£, in two-field inflation models. We use numerical methods to explore the sensitivity of /nl, 
and r to the reheating process, and present simple qualitative arguments to explain our results. 
In general, if a large non-Gaussian signal exists at the start of reheating, it will remain non-zero 
at the end of reheating. Unless all isocurvature modes have completely decayed before the start of 
reheating, we find that the non-linearity parameter, /nl, can be sensitive to the reheating timescale, 
and that this dependence is most appreciable for 'runaway' inflationary potentials that only have 
a minimum in one direction. For potentials with a minimum in both directions, /nl can also be 
sensitive to reheating if a mild hierarchy exists between the decay rates of each field. Within the 
class of models studied, we find that the spectral index ri( , is fairly insensitive to large changes in the 
field decay rates, indicating that uq is a more robust inflationary observable, unlike the non-linearity 
parameter /nl- Our results imply that the statistics of (, especially /nl, can only be reliably used to 
discriminate between models of two-field inflation if the physics of reheating are properly accounted 
for. 
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I. INTRODUCTION 

Inflation has become the dominant paradigm for explaining the generation of the primordial density perturbation 
C, that seeded structure formation, and the Cosmic Microwave Background (CMB) anisotropies. According to the 
standard inflationary scenario, the universe underwent an early period of supcrluminal expansion, stretching the 
primordial density perturbations that were generated by vacuum fluctuations of one or more light scalar fields, 
beyond the causal horizon. On each scale, these fluctuations were promoted to classical perturbations around the 
time of horizon exit. Over time they were gravitationally amplified, and eventually re-entered the horizon laying the 
foundations of all cosmic structure that we observe in the universe today. 

Extending the initial work of Guth pQ, the simplest inflationary mechanism invokes a single scalar field whose 
associated potential has a region which is sufficiently flat to sustain at least 60 e-folds of accelerated expansion [2H1] , 
required to solve the horizon, flatness, and relic problems (see, e.g., [H[U[S]). Whilst single-field slow-roll inflation 
models are consistent with current observational data, there are many reasons to believe that inflation could have 
been driven by more than one scalar field: theories beyond the standard model of particle physics such as string 
theory, supergravity and supersymmetry, generically contain multiple scalar fields. Furthermore, with the possibility 
of greatly enriched field dynamics, multi-field models can give predictions for key physical observables that may be 
quite different from single field inflation models, and thus offer the chance of being constrained. 

Over the last decade, non-Gaussianity has emerged as a powerful probe that may be used to discriminate between 
different models of inflation. Once the power spectrum of £ is known, the assumption that the perturbations are 
Gaussian makes it possible to specify all the properties of the distribution. Any information contained in the depar- 
ture from a perfect Gaussian, non-Gaussianity, is encoded in higher-order correlation functions. Any detection of 
primordial non-Gaussianity, quantified using the non-linearity parameter /nl, would rule out the simplest models of 
single field inflation. 
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A plethora of different mechanisms for generating a large /nl have been proposed in the literature. If the inflaton 
field has canonical kinetic terms then its perturbations are almost exactly Gaussian at Hubble exit and so any 
significant non-Gaussianity must be generated on super-Hubble scales [51 H] • Features in the inflaton potential [5] , 
the curvaton scenario [5HT5]. modulated reheating/preheating [H}|21| . and an inhomogeneous end of inflation [2"2l |2"3"] 
all generate a large non-Gaussian signal. It is also possible to generate significant non-Gaussianity during multi field 
inflation [24, 25 , for a review, see |26j . 

Regardless of the inflationary model, or how many scalar fields were present during inflation, the universe must 
eventually evolve to the hot radiation dominated era of the standard Big Bang model. By the time inflation has 
ended the universe is typically in a highly non-thermal state 1 : the superluminal expansion required to homogenise 
the universe effectively leaves the cosmos at zero temperature, and so a consistent theory of inflation must also 
explain how the cosmos was reheated. This process, which involves a transfer of energy from the inflating field(s) to 
the standard model particles, can be very complex and may proceed via a number of different mechanisms depending 
on the inflationary theory. 

One of the dawning realisations over the last two decades has been that the process by which the universe is reheated 
can have a major impact on physical observables, such as the non-linear parameter /nlj the spectral index and the 
tensor-to-scalar ratio r, that are predicted by the preceding inflationary phase 30 -33J . Indeed, a number of recent 
authors [24j I34H40] have cautioned that the physics of any subsequent reheating phase may affect the observational 
predictions of their inflationary models. Generically, we should expect the inflaton to couple to other fields which do 
not play any role in driving inflation, and such interactions are unavoidable from an effective field theory perspective. 
For example, it has been shown that the inclusion of such interactions can lead to particle production effects, which 
radically modify the phenomenology of some inflation models [4ITI44] . 

As emphasised in [3S], to connect the physics of inflation with observation, the statistics of £ should ideally be 
followed all the way up until the time of last scattering, where the microwave background anisotropy was imprinted. 
Without a fundamental UV complete theory describing all early universe physics, this is unfortunately impossible. 
Thankfully however, we may rely on the fact that in the absence of isocurvature (entropy) modes, the curvature 
perturbation becomes a conserved quantity on superhorizon scales. This was demonstrated to all orders in cosmological 
perturbation theory [JS] and was also verified using a gradient expansion method [3H1 HZ] • Hence, the statistics of £ 
evaluated when all isocurvature modes are exhausted and an adiabatic condition reached, are those that are measured 
today. 

Many previous works have assumed that the universe is reheated instantaneously [48H51) . fossilising the curvature 
perturbation immediately. But this is an idealisation: reheating presumably takes a finite time to complete. Recently, 
the authors of [52 have shown explicitly that for canonical single field inflation models with quadratic minima, an epoch 
of preheating does not alter the amplitude of the scalar bi-spectrum generated during inflation. Ideally, the evolution 
of ( should also be followed through the subsequent phase of reheating where the energy of the oscillating inflaton 
is transferred to radiation. However recent studies have shown that the amplitude of the curvature perturbation 
remains unaffected even during perturbative reheating, see for example |53j . Thus naively, one might expect the 
scalar bi-spectrum to also remain unchanged by this process since £ itself is conserved at a non-linear level for single 
field models, although this remains to be seen explicitly. 

Even less clear is how reheating affects the evolution of C in multi-field models, since when more than one field is 
present, isocurvature fluctuations can cause C to evolve on super-Hubble scales. It is already known that the two- 
point correlation function of £ can be affected by metric preheating |54j . Until an adiabatic condition is reached, such 
as in the case when the universe is radiation dominated, all observable quantities associated with £ continue to evolve. 
How sensitive then are the key inflationary observables to the reheating process? Is this sensitivity heavily dependent 
on the inflationary model? Does the level of non-gaussianity that exists at the end of inflation survive until the 
completion of reheating? Focusing on two-field inflation models and assuming that reheating proceeds perturbatively, 
it is the purpose of this paper to address these questions. 

By numerically implementing the SN formalism, we follow the evolution of £ beyond the end of inflation, until the 
completion of a phase of perturbative reheating. We parametrise the decay of the oscillating inflation and isocurvature 
fields into relativistic particles by introducing decay terms into the field equations. Our goal is to investigate the 
sensitivity of the key inflationary observables to the physics of perturbative reheating. We study two classes of 
potential: the 'runaway' type which has a minimum in only one direction; and potentials which have a minimum in 
both directions. 

The paper is organised as follows: in Section ([TTJ) we recall the SN formalism and review the textbook elemen- 
tary theory of reheating, before discussing its numerical implementation within the separate universe picture. In 
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Section (III I we study the evolution of /ml and other £-related statistics during the reheating phase for the class 
of potentials which posses a minimum in only one direction. Then, in Section (IV I we repeat the same analysis for 
potentials where both directions have a minimum. We also consider an example of non-separable potential models 
in Section ([V]). We discuss and conclude in Section (VI I. The expert reader familiar with the elementary theory of 
reheating and the SN formalism may wish to omit Sections (II A I and (II B). 



II. PERTURB ATIVE REHEATING, NON— GAUSSIANITY AND THE <5N FORMALISM 

The two-field inflation models that we study in this paper are described by the action 



S = I d x^f^g 



(1) 



where M p = l/y/8nG is the reduced Planck mass. The standard slow-roll parameters are defined as 
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where the subscripts denote differentiation with respect to the fields. 

We will consider various forms of W(ip, x)i with the only constraint that W(tp, x) must have a minimum in one, or 
both of the field directions to enable at least one field to oscillate and reheat the universe. The fundamental difference 
in form between one-minimum and two-minima potentials provides a logical division of our analysis into classes. This 
is partly motivated by the work of [SS] where two broad classes of behaviour for the evolution of Q were recognised: 
Potentials that contain a 'natural focussing region', which is guaranteed for a two-field model with minima in both 
directions, allow neighbouring trajectories in field space to converge 'naturally', quenching the flow of power from 
isocurvature modes to C- Alternatively no such focussing region may exist, which is the case for a two-field model 
with only a single minimum, and so Q will continue to evolve until an adiabatic condition is reached. In the latter case, 
predictions for observables such as /nl cannot currently be linked directly to the physics of the inflationary model, 
as they will be dependent on the subsequent phase of reheating. Even in the former case, if the universe approaches 
adiabaticity by the inflating/isocurvature trajectories converging in, and oscillating about, their global minima, then 
it is not clear how the decay of the oscillating fields into radiation affects the final stages of the evolution of £■ We 
note that adiabaticity may also be reached via a third waterfall field, as is the case in hybrid inflation |55) . 

In the following subsections we introduce the simple perturbative reheating scheme that we use throughout this 
paper, and briefly review the <5N formalism that is used to compute the statistics of £. 



A. Elementary theory of reheating 

In this section we recall the elementary theory of reheating based on perturbation theory that was developed 
in |56[ 157] . For the purposes of this discussion, we assume for simplicity that the potential W(ip,x) m the action 
Eq. ([l]) only has a single minimum in, say, the x direction, and we assume that this is the inflationary direction. 
The reheating mechanism presented here only applies to the directions in the potential that are associated with 
well-defined minima. At inflationary energy scales we may neglect the contribution to gravity from any other fields 
such as Xb bosons (not to be confused with the inflationary x field) or ifjf fermions in the action Eq. (JlJ. Hence, for 
cosmological applications we may retain only the dominant fields ip, x and gravity. Then in a flat FRW universe, the 
Friedmann equation reads: 



H 2 



I 



3Af p 2 



+ \x 2 + W{if, X ) 



(3) 



The dynamics of x is governed by the Klein-Gordon equation 



X + ZH X + W, X =Q, 



(4) 



4 



and similarly for tp. For sufficiently large initial values of XiV > -^pj Hubble friction dominates over x (and tp) and 
the potential term W((p, x) in Eq. ^ is assumed to dominate over the kinetic terms. During this slow-roll stage, 
the universe inflates, expanding quasi-exponentially. As the inflating x held rolls toward its minimum at Xo it gains 
kinetic energy, eventually bringing inflation to an end, whilst the tp field continues to contribute to the expansion rate. 
We assume that the minimum in the x direction is quadratic to leading order, ^m x x 2 , and so xo — 0. We note that 
a similar discussion may also be applied for theories with quartic minima |58j . Ignoring for the moment the effects of 
particle production, as the inflaton approaches and inevitably overshoots its minimum, it begins to oscillate about xo 
on a shorter time scale compared to the Hubble time. Here, we assume that H <C m x after inflation has ended. 2 The 
frequency of the oscillations is ko — m x . The large vacuum energy of the inflaton then exists in spatially coherent 
oscillations, which can be interpreted as a collection of a number of x _ particles with zero momenta. The density, 
n x = p x /m x of this coherent wave of particles decreases as a~ 3 , since the condensate behaves as non-relativistic 
pressureless matter: p x — \{x 2 + m x X 2 ) ~ a ~ 3 - 

The amplitude of the x oscillations gradually decays due to the Hubble expansion and also because of the transfer 
of energy to lighter particles produced by the oscillating field. As these decay products thermalise, the Universe is 
reheated. The inflaton may decay into bosons Xb and fermions tbt due to —\g 2 X 2 x\ and —htpfipfX interaction terms, 
which should now be included into the fundamental action Eq. (fTl). Based on the above interpretation of the spatially 
homogeneous, coherently oscillating x field, the effects of particle production may be incorporated into Eq. Q |59j : 

X + 3H X +{m 2 x +n(k ))x = 0- (5) 

Here, n(/c ) is the flat space polarisation operator for the field x at four-momentum k = (fc , 0, 0, 0). It can be shown 
that the real part of Tl(ko) gives only a small correction to m 2 , but when fco > min(2m Xi , 2m,/,,), n(fco) acquires 
an imaginary part Imll. Working in the limit m x 3> {H, Imll}, which are conditions that should be satisfied after 
inflation, neglecting the time-dependence of Imll and assuming H = 2/3i, the approximate solution to Eq. ^ is: 

X (t) « Jf p exp (~Tt) sin (m x t) , (6) 



37rm x f \ 2 

where Y = T(x XbXb) + r(x — > ipfipf) is the total decay rate of x particles. Here we have used the relation 
Imll = m x T which follows from unitarity [3D]- Eq. ([6| implies that the amplitude of the \ oscillations decays as 
X (t) ~ a- 3 / 2 exp(-irt). 

For a phenomenological description of the reheating effect, one can add an extra friction term T x x to the classical 
equation of motion of the field x, instead of adding the polarization operator j581 161] : 

X+(3H + T x )x + W tX = 0. (7) 

Once again assuming H = 2/3t and a quadratic minimum, W — ^m 2 x 2 , the solution of this equation is exactly 
Eq. (p3l). Multiplying through by x it is intuitive to rewrite Eq. (T7J) as p x + 3Hx 2 + Y x x 2 — 0. Now, since x is rapidly 
oscillating around xo approximately sinusoidally, it can be replaced by its average over a single oscillation cycle 3 , 
(x 2 ) C ycic = Px- If the decay products of the oscillating x field are very light relative to x itself, and are only bosonic, 
we can model them as a (single) rclativistic radiation fluid: 

p 7 + \H Pl = Y xPx = T x x 2 , (8) 

H 2 = (P X +P V + Ar) • ( 9 ) 

At this point a number of comments surrounding the validity of Eqs. ([7| and ^ are in order. Firstly and most 
importantly, this simple phenomenological equation ([7]) is only valid when x is rapidly oscillating about Xo : the 
'particle creation' term, T x x, should not be present during inflation. Furthermore, since in this example the tp field 
does not have a minimum about which it can oscillate, it should not be coupled to radiation: T v = 0. Secondly, Eq. ^ 
(as is Eq. fcfy) is valid only when m x 3> H and m x >• T x . We have also made the assumption that the decay rate T x of 



2 For potentials with local curvature much different to ijn^x 2 , this estimate can be very different. 

3 If the motion of \ is approximately that of a simple harmonic oscillator, (V) = (x 2 /2) = p x /2 an d so we see that (P x ) = (x 2 /2 — V) 
vanishes and the coherent oscillating \ behaves as pressureless matter, justifying our previous statements. 
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the inflaton can be calculated using the standard methods of quantum field theory, describing the decay x — ► XbXb- If 
however, many Xb _ particles were produced in the early stages of particle production, the probability of decay becomes 
greatly enhanced by effects related to Bose-statistics, which may lead to explosive particle production [52"HS4"] . When 
the amplitude of the oscillating field is sufficiently large, we should also expect reheating to occur in a different way 
through parametric or stochastic resonance |64H66j . 

In this perturbative scheme, reheating completes at time t c , when the Hubble rate H 2 = p/3M 2 ~ t~ 2 drops below 
the decay rate T x . The density of the universe at this moment is then 

p(t c ) ~ 3H 2 (t c )M 2 = 3T 2 X M 2 . (10) 

If the decay products interact with each other strongly enough, then thermal equilibrium is quickly established and 
may be maintained at a temperature Tr. Treating this ultrarelativistic gas of particles with Bose-Einstein statistics, 
the energy density of the universe in thermal equilibrium is then 

p(T R )~(^)g.l£, (11) 

where the factor <7*(Tr) ~ 10 2 — 10 3 depends on the number of ultrarelativistic degrees of freedom. Comparing 
Eqs. (10) and (11) we arrive at 



TR-OTv/iyWp. (12) 

In order not to spoil the success of BBN, the inflaton decay products should be quickly thcrmalizcd through scatterings, 
annihilations, pair creation and further decays, such that the universe is completely radiation dominated before 
the BBN epoch. This constrains the reheating temperature to be Tr > 5MeV [571 HEj, which in turn implies 
r x > 4 x 10 _40 M p . We ensure that this bound is always satisfied throughout our paper. For such weak decay rates, 
reheating would proceed incredibly slowly if the process were entirely perturbative. In reality however, as alluded 
to above, the universe is unlikely to be reheated via a mechanism that can be described completely by standard 
perturbation theory, and so we interpret such bounds on T x rather loosely. There is also an upper bound on Tr (and 
so T x ) coming from the overproduction of gravitinos [69| 170] , which does not apply to us as we are not considering 
supersymmetric models. 

Despite various limitations, the elementary theory of reheating is appealing due to its simplicity and ability to be 
very successful in describing the reheating process in certain regimes. In this paper we are interested in the effects that 
reheating has on the evolution of statistical properties of £, such as /nl- To this end, we parametrise the reheating 
process with Eqs. Q and ^ and assume that any important physics that may affect the evolution of £ are well 
described by this paramctrisation. Since we are concerned with two-field models of inflation, we also assume that 
this description of reheating applies to both fields, ip and %, subject to the limitations discussed above. 

Whilst reheating may well be more complex than the simple perturbative model we consider, it is a useful scheme 
for determining how sensitive the primordial observables may be to reheating, and to check whether any general trend 
exists across different models. For example, one might speculate that any large non-Gaussianity is generically damped 
to zero by reheating, as is often (but not always |71j ) the case during inflation if the isocurvature mode decays during 
slow-roll [36j[72]. We will show that this is not the case for reheating. 



B. The <5N Formalism and non— Gaussianity 

The SN formalism [73H75] has been used extensively throughout the literature to compute the primordial curvature 
perturbation and its statistics. The formalism relates £ to the number of e-folds of expansion N, given by: 

N(U,t c ) = J° H(t)dt, (13) 

which is evaluated from an initial flat hypersurface to a final uniform density hypersurface. The perturbation in the 
number of e-foldings, 8 N, is the difference between the curvature perturbations on the initial and final hypersurfaces. 
We take the initial time to be Hubble exit during inflation, denoted by i*, and the final time, denoted by t c , to be a 
time deep in the radiation dominated era when reheating has completed. The curvature perturbation is then given 
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by [75] (or [7B] for the covariant approach) 

C = SN = NjSipj, + N ,u 5( Pi* 5 VJ* + ■ ■ ■ , (14) 
/ ij 

where N,j = dN/(dipl) and the index I runs over all of the fields. In general, N(t c ,t*) depends on the fields, <fii(t), 
and their time derivatives, (fi{t). However, if the slow-roll conditions, 3Hipj ~ —Wj, are satisfied at Hubble exit, 
then N depends only on the initial field values. The radiation fluid remains effectively unperturbed at horizon exit 
as it does not yet exist, and so does not feature in the above expansion. The power spectrum and bispectrum defined 
(in Fourier space) are given by [5]: 

2tt 2 

(C kl Ck 2 ) = (2^) 3 ,5 3 (k 1 +k 2 )^ c (fc 1 ), (15) 

(CkiCk.Cka) = (27r) 3 S 3 (k 1 + k 2 +k 3 )B ( (k ll k 2 ,k 3 ). (16) 

From this we can define three quantities of key observational interest, respectively the spectral index, the tensor-to- 
scalar ratio and the non-linearity parameter 

»t-l = in— 77' ( 17 ) 



/NL 



dlogk 

5 fc 3 fc 3 fc 3 B c (fci,fc 2 ,fc 3 ) 

6 k\ + fc 3 + fc| 4tt 4 ^ 



(18) 
(19) 



Here V* is the power spectrum of the scalar field fluctuations and Vt = 8"P* = 8ff^/(47r 2 Mp) is the power spectrum 
of the tensor fluctuations. As defined above, /nl is shape dependent, but it has been shown that the shape dependent 
part is much less than unity [6J [3 [77H79] for local non-Gaussianity in canonical models. Since the ideal CMB 
experiment is only expected to reach a precision of /nl around unity [80] . we calculate the shape independent part 
of /nl, denoted by / N ^ in [ZZ1 [HI] - Whenever the non-Gaussianity is large, |/nl| > 1, as is the case considered 
throughout this paper, we can associate / N ^ — /nl- This k independent part of /nl and the spectral index can be 
calculated by the SN formalism, 

v c = E^*' ( 2 °) 
i 

nc-l = -2e* +7 f -^ 'V./' ■' . (21, 



2 Eij<P*jNj iM,i 
H* EkN 2 k 



/nl - g -2 — ■ (22) 

We use the same sign convention for /nl as the WMAP team [82 . The latest observations from 7 years of WMAP 
data are [55] 

n c = 0.967tg;g^ (assuming r = 0) , (23) 
r < 0.36 (95% CL) , (24) 
-10 < /^jf 1 < 74 (95% CL) . (25) 

The crucial difference between single and multi-field inflation is that in single field inflation, the slow-roll solution 
forms a one-dimensional phase space. Hence, by virtue of the attractor theorem there is a unique inflationary 
trajectory that is always quickly reached. Furthermore, the end of inflation takes place at a fixed value of the 
inflaton field, corresponding to a fixed energy density. When two fields are present however, the phase-space is two- 
dimensional with an infinite number of possible classical trajectories in field space. The values of the two fields at the 
end of inflation will in general depend on the choice of trajectory. Then, to compute the SN derivatives (Nj etc) in 
multi-field models, an extra piece of information, a conserved quantity along a given trajectory, is required. Within 
slow-roll, such a constant of motion exists (see for example [771 US] ) > and under the assumption that the potential 
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is sum separable W = X(U(p) + V(x)) a , or product separable W — \(U (ip)V (x)) a in the fields, explicit expressions 
for the SN derivatives may be obtained p5 l l36 l [77 | I5T] . 

Recently, the authors of [53] used raytracing techniques to reformulate inflationary perturbation theory in the 
language of geometrical optics. Whilst this technique yields differential equations from which the SN coefficients 
can be computed efficiently, closed-form expressions for the resultant path-ordered exponential integrals can only be 
obtained under conditions of separability and slow-roll. By decomposing the field perturbations into curvature and 
isocurvature perturbations, similar expressions are also found in |85j . 

Exact solutions, valid beyond slow-roll, have been obtained assuming a sum-separable ansatz for the Hubble 
parameter [86 . Such an ansatz, besides being very restrictive, cannot be applied to a phase of perturbative reheating 
as it relies upon monotonicity of the field variables. 

Alternative long-wavelength (LWL) formulae have also been developed to analytically study the nonlinear evolution 
of long wavelength cosmological perturbations in the early universe |87j . In such an approach, the perturbations are 
written in terms of quantities of the corresponding exactly homogenous universe to the leading order of the gradient 
expansion. The formulae have recently been extended to study nonlinear perturbations in universe where multiple 
scalar fields and perfect fluids coexist 88, 89 . 

In this paper, we will use the SN formalism. In all cases, what currently evades us is a practical method for 
analytically computing the SN derivatives for arbitrary potentials during slow-roll and separable (and non-separable) 
potentials beyond slow-roll. While still assuming slow-roll at Hubble exit however, one can go beyond the slow-roll 
approximation and the condition of separability by numerically solving the second order equations of motion, Eqs. ^ 
and (TSJ) , together with the Friedman constraint J9J, introducing the decay terms T v and T x when applicable. 



C. Numerical Code 



The SN formalism is based on the assumption that (smoothed) spatially separated patches of the universe will 
evolve on super-horizon scales like independent, unperturbed universes up to small corrections. This is the separate 
universe picture [SUl EE]- An ensemble of smoothed regions picks out a collection of trajectories in phase space 
which is often referred to as a 'bundle' [38, 84J. In essence, the SN formalism requires knowledge about how such a 
bundle, centred on a fiducial trajectory, evolves. Our choice of gauge demands that each trajectory in the bundle is 
evolved from an initially flat hypersurface up to a hypersurface of constant energy density. Hence, each trajectory will 
experience a slightly different expansion history in order to bring them to a common energy density. The adiabatic 
mode is generated by fluctuations along the fiducial trajectory, whilst fluctuations between neighbouring trajectories 
generate the isocurvature modes. 

Acknowledging this simple picture, the SN formalism may be implemented numerically: First, the fiducial trajectory 
emanating from {p*,x*} is constructed by solving the full, non-linear second order field equations, i.e., Eq. Q (and 
similarly for the p field) together with the Friedmann constraint ^ and the equation for the radiation fluid (|8j) . The 
bundle is then formed by evolving neighbouring trajectories with slightly perturbed initial conditions, ip* — > p* + Sip* 
and x* ~ * X* + $X*- Each trajectory in the bundle is then brought to a common energy hypersurface where the 
partial derivative of N(t c , t*) with respect to the field values at horizon crossing {ip* , \*} is taken using a seven-point 
'stencil' finite difference method |92j . This provides a fast, efficient method for computing r and /nl for an 
arbitrary two-field model, valid beyond slow-roll and through a phase of reheating. Numerical codes based on the 
moment transport equations have also been developed |93j . 

As discussed in Section (II A), the reheating parameters T v and r x are set to zero during inflation. It is only when 
each individual trajectory in the bundle passes through its minimum {xo, <Aj} for the first time that T v and T x are 
introduced to the field equations, sourcing the radiation fluid. In general, for any given trajectory, p will not reach 
the minimum of its potential at the same time as x, and so and T x are 'switched on' at different times along the 
same trajectory. Furthermore, for each directions in the potential, the foliation of the entire bundle of trajectories as 
determined by each trajectory reaching xo (and likewise po) does not in general occur at a surface of constant time 
or a surface of constant energy, but rather at a surface of constant xo (and <po) 4 - We refer to these surfaces as the 
reheating hyper surf aces. For potentials which have minima in both directions there are two such hypersurfaccs. If the 
potential does not have a minimum in the x (° r f) direction, then T x = (or T v = 0) always. Furthermore, we also 
ensure that when the potential has a minimum in, say, the x direction, the conditions m x 3> T x and m x 3> H are 
satisfied. This definition of the reheating hypersurface is more refined than that of [38], where reheating was initiated 



4 This is true for global minima. If the oscillations of one field, \ sa yi occurred in a local minimum, which is a function of the other field, 
Xo(<p), this statement will not hold true. We do not consider such models in this paper. 
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at a surface of constant density. It is also different to that of [S3], where the decay terms were present throughout 
inflation. The rest of this paper is dedicated to exploring the sensitivity of /nl, n( and r to the reheating process. 



III. ONE MINIMUM 



In this section we present numerical results for the statistics of £ for the class of two-field potentials which have 
one minimum in the x direction. In what follows, \ may be identified as the inflaton and p as the field which sources 
the isocurvature perturbations. The ip field is not directly involved in the reheating phase and so T v = at all times. 



A. Quadratic minimum: W(ip, x) = W / oX 2 e Xv 

This potential was first introduced by [55], and has made frequent appearances in the literature since then [58"1IT2"II55T - 
[55] , It does not contain a 'focussing' region where neighbouring trajectories in the bundle may converge, due to its 
'runaway' form in the ip direction. Hence, C an d its statistics will continue to evolve after inflation has ended. 
The parameter space for which /nl may be large at the end of inflation was derived in [25]. Essentially, the initial 
background trajectory must be fined-tuned to be nearly parallel to the axis of the inflaton. It is useful to first consider 




r x =o 



55 57.5 60 62.5 65 67.5 70 64 66 68 70 72 

N N 



FIG. 1: Potential: W{ip,x) = Wa^e~ Xifi ■ Left panel: The evolution of the background fields for A = 0.05, ip* = 10" 3 Af p and 
X* = 16.0M P . Right panel: /nl as a function of TV (with T x — 0) after inflation has ended. The parameters used are: A = 0.06, 
(p t — 10 _3 M P and x* ~ 16.0M P . In both panels, the solid vertical (black) line denotes the end of inflation, 7V C . 



the evolution of the fiducial background fields and /nl in the limit of no reheating, i.e., T x = 0. We set A = 0.06, 
<p* = I0~ 3 M P and x* = 16.0M p . With this choice of parameters, a large /nl is still present as slow-roll breaks down 
(at N ~ 64.5). However, since no limiting trajectory is available, /nl continues to evolve: as x reaches the minimum 
of the potential, it undergoes damped oscillations about x = 0, which induces large oscillations in /nl- The ip field 
remains in slow-roll after inflation has ended and throughout the period of oscillations of x an d continues to evolve 
towards ever increasing values. We note that if the potential is sufficiently steep in the tp direction, i.e., large A, then 
the slow-roll conditions for ip may be violated. Such large values of A tilt the initial trajectory away from the axis of 
inflaton and hence a large /nl cannot be generated. As such, we do not consider such regions of parameter space. In 
Fig. [T] we show the late time evolution of the scalar fields and /nl- The 'spikes' in the oscillations of /nl correspond 
to the x field changing direction at the maximum of its oscillation. Fluctuations between neighbouring trajectories in 
the bundle source /nl- These trajectories continue to diverge in the tp direction due to the geometry of the potential 
in this region and so /nl continues to evolve, decaying towards zero. If the evolution of /nl were followed indefinitely 
with T x = 0, we should expect it to settle to f&if 1 ~ 0. In order to explain this, we need to examine each derivative 



term contributing to the expression for /nl , Eq. ( 22 ) . We begin by considering the slow-roll solution for tp and the 
r) vv slow-roll parameter, 

ip = p,e 2XN , Vvv = 2A [2\<ple 4XN - l] , (26) 
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0. This slow-roll solution 



which shows that trajectories will continue to evolve indefinitely in the tp direction if T x 
for tp is shown against the exact numerical solution in the left panel of Fig. [TJ 

Now, inspecting the individual derivative terms (iVj = ON/ (dpi) etc) contributing to /nl 
practically constant, N x ss (2e*) -1 / 2 , throughout the entire inflationary and post-inflationary phase, acquiring this 
value when the fields leave the horizon. To explain this, we assume that H is monotonic in time, enabling us to 



wc find that N x remains 



re- write Eq. ( 13 ) as 



N ■ 



dH 2 
2H 



(27) 



Taking the derivative with respect to x* we find 

, 1 OH 2 

A v 



2H dx* 



d 

dx* 



2H 



AH 2 



(28) 



where the derivative inside the integral is computed by holding H constant. The derivative at the boundary c vanishes, 
since by definition the surface c corresponds to one of constant H. Using the fact that the fields are in slow-roll at 
horizon exit, the first term on the RHS of Eq. (28 1 reduces to (2e* )~ 1 / 2 . Then, to explain why N x remains constant at 
this value requires arguing that the integral term m Eq. ( 28 ) is negligible, i.e., after perturbing x* > surfaces of constant 
H must coincide with surfaces of constant H. This is indeed the case if a hierarchy of kinetic energies exists between 
the fields at horizon crossing, i.e., |x*| 3> \p*\- Since the kinetic terms arc canonical, the fields follow the gradient of 
the potential, and as they are in slow-roll at horizon exit, this hierarchy implies |W X |* 3> (W^l*. If this is the case, 
the dependence of H on x* is rapidly washed out, and the two-dimensional bundle in the x direction (holding tp* 
fixed) degenerates to a caustic. We have found that the condition \W X \* ^> |W V |* is sufficient to guarantee that the 



integrand of Eq. (28 1 is always small from horizon crossing until oscillations of \ begin. During the oscillatory phase, 



the integrand oscillates about zero with an amplitude that decays with the Hubble expansion, and when integrated 
over many oscillations, the net result is a negligible correction to N x . By the same argument, N xx remains roughly 
constant at N xx w 1 — (r] xx /2e x )^, which, for this particular potential is independent of A and the field values at 
horizon crossing, N xx |. 

Furthermore, we find that the following approximate scaling relations hold to a remarkable accuracy throughout 
the entire inflationary and post-inflationary evolution: 



N ir , 




(29) 
(30) 



The scaling relation between N vip and N v was first derived in 



a 'ridge', situated at tp = 0, of a generic potential. 



by considering a first order Taylor expansion about 
Assuming the slow-roll conditions, the same analysis applies to the 



model we study here as long as the potential remains well approximated by W ~ WqX (1 — )j i- e -> higher order 



terms in \tp remain small. This requires « C(A' 
of trajectories rolls off the ridge: tp — tp*e a ( H *~ H \ 



, In this regime, tp grows exponentially with H as the bundle 
3 A/2 Wq . A short calculation reveals 



(31) 



where j3 is some model-dependent constant. We refer the re ade r to |38) w here the complete derivation is presented 



Taking g^- (on surfaces of constant H) on both sides of Eq. (31 ) gives Eq. (29). Similarly, taking the derivative with 



respect to x* an d using the tp slow-roll solution Eq. ( 26 ) gives Eq. ( 30 ) . 

We show evolution of the N v , N vv and N vx derivatives before and after inflation in Fig. [51 which clearly illustrates 
the scaling behaviour captured in Eqs. (29) and (30 1. Remarkably, not only does this scaling behaviour hold after 
inflation has ended, but it also holds during reheating. Indeed, we find that it remains an excellent approximation 
across the entire range of T x that are within our numerical capabilities, including T x = 0. 

The derivation of these scaling relations as sketched above relies on a number of approximations, including slow- 
roll. The sub-dominant field tp always remains slowly rolling, however \ necessarily does not. This does not seem 
to violate Eqs. (29) and (30), suggesting that validity of these relations are more reliant on tp being a linear function 
of tp* , and that tp grows exponentially as the bundle slides off the ridge. As mentioned above, these conditions will 
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FIG. 2: Potential 
Evolution of the derivatives N vv 



W{tp, X ) = W 0X 2 e 



Numerical verification of the scaling relations Eqs. (291 and (301 



and N v . The horizontal dashed line in the lower panel represents the value of tp* 



Left panel: 
the constant 



of proportionality between N vv and N v . Right panel: Evolution of the derivatives N vx and N v . The horizontal dashed line in 
the lower panel represents the value jx(W x /W)* = j x (2e x ) 1 ^ 2 , the constant of proportionality between N vx and N v . We show 
evolution of the derivatives for the last few e-folds of inflation, up until ( has become conserved at the completion of reheating. 
We see small departures from scaling at the start of reheating as \ oscillates about its minimum, but as \ settles down, the 
scaling behaviour is quickly recovered. In both panels, the parameters used are: A = 0.05, to, = 10~ 3 M D , y* = 16.0A/ D and 
F x = \/10 _1 WoMp. The solid vertical (black) line denotes the end of inflation, N e 
the start of reheating, N T . The Hubble rate at the start of reheating is H r « y/7 x 10~ 2 WqM p . 



tp, = ru "M p , x* - ^.^., p 
and the dashed vertical (blue) line denotes 



break down when tp ~ 0{\ 1 I 2 ). Then, using tp ~ A */ 2 in Eq. (26) we may very roughly estimate how many e-folds 
we expect the scaling relations to remain valid: N ~ ^ln (A -1 /^/?*). For example, for A = 0.05 and tp* = 10~ 3 we 
have N ~ 85. 



We now return to the expression for /nl, Eq. (22). The N xx and N vx derivatives are in fact negligible compared 
to N vtp and so can be safely neglected. Making use of the approximations discussed above, we may write /nl solely 



in terms of N v : 



/NL 



N 3 

(32) 



where = N x « (2e x ) 1 ^ 2 . The asymptotic behaviour of /nl is clear: trajectories in the bundle continue to diverge 
away from one another in the tp direction according to Eq. (26), which continuously sources N v , making it grow 



increasingly more negative. Hence, in the limit that N v — > — oo we expect /nl — ^ 0' which is what is observed in 
the right panel of Fig. [TJ justifying our previous claims. The sign of N v can be argued from the geometry of the 
potential: diverging trajectories source negative N v (as the sign of Eq. (31 1 indicates), whilst converging trajectories 
source positive N v 38 . 

With the limiting case r x = understood, we now move on to explore the dependence of /nl, 3,1 on T x , keeping 
the same parameter choice A = 0.06, tp* = 10 _3 M p and \* = 16.0M p . We express the decay rates in terms of the 
overall potential normalisation, Wo . Whilst its value sets the scale of inflation and determines the amplitude of the 
primordial power spectrum and hence is constrained, it does not affect the statistics of C and so we leave Wo as a free 
parameter. Where applicable, we also give the value of the Hubble rate at the start of reheating, H r , in units of Wo 
so a direct comparison between the expansion and decay rate can be made. 

The switching on of the decay terms at the reheating surface sources the radiation density. As the \ field oscillates 
about its minimum, its kinetic energy is transferred to the radiation fluid, resulting in bursts of particle production. 
As radiation fills the universe, Hubble damping slows the motion of tp to a crawl and as we approach fi~ ~ 1, it 
asymptotes to a constant: tp(t — > oo) ~ const. Herein is the fundamental difference in the motion of tp when T x 7^ 
compared to Y x — 0: as radiation comes to dominate, trajectories in the bundle cease to evolve. The bundle does not 
degenerate to a caustic as would be the case if the trajectories were naturally focussed by a region of the potential, 
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FIG. 3: Potential: W(ip, x) = W / oX 2 e _A¥> • We show the evolution of /nl during reheating for various decay rates F x , which are 
in units of vWoM p . In both panels, the solid vertical (black) line denotes the end of inflation, N c , and the dashed vertical (blue) 
line denotes the start of reheating, N r . Left Panel: The parameters used are: A = 0.06, ip, = 10" 3 M P and x* = 16.0M P . The 
Hubble rate at the start of reheating is H r « V7 x 10~ 2 WoM p . Right Panel: The parameters used are: A = 0.05, if, = 10 _3 M P 
and x* = 16.0Mp. The Hubble rate at the start of reheating is H r « y/6 x 10- 2 VK A/ P . 



but nonetheless this freezing of the ip field guarantees that £ becomes conserved. This does not happen in the T x = 
limit where the trajectories continue to diverge in the tp direction, always sourcing £. 

In the left panel of Fig. [3] we show the final stages in the evolution of /nl as a function of iV for various decay rates 
T x . Most importantly, we see that reheating does not damp out /nl to zero. We interpret the fine details of the 
plot as follows: At the end of inflation (N c = 64.56) a large, negative /nl is still present, and just before reheating 
begins 5 (N T — 65.10) /nl is growing increasingly more negative. We see that as the decay rate F x is increased from 
zero, I/nl^I freezes out to larger values. In another example where /nl is decaying toward zero as reheating begins, 
the effect of increasing the decay rate from zero is to freeze out I/nl^I to smaller values. This is shown in the right 



panel of Fig. [3] 

.1 ,,,! ,,r l , 

/NL 

of /nl on N v . Let us begin by considering the splitting 



This opposite dependence of I/nl^I 011 T x for A = 0.05 and A = 0.06 is a consequence of the non-trivial dependence 



J* 2H J* 2H J r 2H 

Here Nq is the number of e-foldings from horizon crossing (t*) up to the start of reheating (t r ) and Ni is the number 
of e-foldings from the start of reheating up to radiation domination (t c ). Firstly, it is important to appreciate that 
A^o contains contributions not only from the slow-roll inflationary phase, but also from the non-negligible post- 
inflation/pre-reheating evolution, that must be accounted for. Whilst the standard methods (see eg. Refs. [25 ) l36l l77l 
[ST]) may be used to compute the derivatives (N j etc) of the slow-roll contribution to Nq, derivatives of the remaining 
non-slow-roll contribution to Nq cannot be calculated explicitly. Secondly, Nq does not contain any dependence on 
the reheating process. Since we are interested here in studying the effects of reheating on /nl 11 , we compute Nq and 
its derivatives numerically and focus on trying to understand the correction Nx, which contains all the dependence 
on T x . 

For the derivative of the correction Ni with respect to (p* we need only consider the term 

N llV — f C 77— (X) dff 2 , (34) 
Jr d<p* \2HJh 



5 Here, we use the terminology 'start of reheating' to refer to the time when the fiducial background trajectory, emanating from {x*> V*}, 
crosses \0 f° r the first time. 
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FIG. 4: Potential: WUp,x) = W x 2 e~ X,p ■ Left Panel: The algebraic function f^ 1 



correction Ni ilp 
right panel. Right Panel 



as a function of the final value of the 
0.05 curve which correspond to the decay rates given in the 
The evolution of the derivative JV V = No.^, + iVi,^ for the same decay rates as Fig. [3] for A = 0.05. 

and the dashed vertical 



Eq. (351. We label the positions along the A 



All decay rates are in units of vWoM p . The solid vertical (black) line denotes the end of inflation, N t 
(blue) line denotes the start of reheating, N T . The Hubble rate at the start of reheating is H r « \/7 x 10~ 2 WoM p 



since the derivative at the boundary at r cancels with the N contribution and derivative at the boundary c vanishes 
since c is defined as a surface of constant H. Since H is a function of l f{t) and Pj(t), all of which depend on 
(/j*, this integral cannot be performed analytically beyond slow-roll. However, we can make progress by using our 
results, N x (2e*) -1 / 2 , {N vx , N xx } « N vv and N vv ~ N^/ip* which also hold during reheating. Then, using 
the fact that during reheating N± tX ~ 0, and taking the time t c to be deep in the radiation dominated era such that 
N\ iV — const, Eq. (32) becomes 



.final _ _5_ (AW + AW) 3 , , 

~ 6\<P*\ [(AW + AW) 2 +3 2 ] 2 ■ {M> 

We plot this algebraic function, f^ff 1 against Ni tV , in the left panel of Fig.jijfor three different choices of the potential 
parameter A — {0.05,0.06,0.07}, with the same field values at horizon crossing (/j* = 10 _3 M p and x* = 16.0M P . 
Changing A obviously changes g* and modifies the evolution of the bundle, changing Nq iV> . I n the right panel of Fig. [4] 
we show the evolution of A^ for various decay rates with A = 0.05. The final values of AAi i¥ , (final) = N v (final) — A^W 
are marked on the corresponding curve in the left panel of Fig. [4] Only the Ni jip < region of Eq. (351 is physical: 
we have already argued that diverging trajectories can only generate negative Ni jtp , and we have confirmed this 
numerically. As can be seen from the left panel of Fig.|4j Eq. ( [35] ) has three stationary points at finite Ni ttp : 

- N 0>v , -AW ± V3g* . (36) 

The N% v — —AW root is an inflection point where /nl 3,1 — 0- The N\ tV = — A<W + VSg* root is a local maximum 
where /nl 3,1 would be always positive and so is not physical. The minimum at JVj„ = —AW — v3<7* however is 



physical and bounds the maximum value of I/nl^'I when Eq. (35) has a minimum at negative A 7 ] 



1 / 75 

|maxR< k^JV 1024' f ° r N ^ + ^9*>0- (37) 

A minimum at negative is clearly seen in the left panel of Fig.|4]for A = 0.05. If on the other hand, the minimum 
exists at positive AW; (i.e., AW + V3g* < 0) then the maximum value of |/NL, al | is instead bounded by its value at 
the start of reheating: 

c TV 2 

|/£rW « \hL(tr)\ w R , | m2 °'% a , for N „ + y/3g. <0. (38) 
"Iv 3 * I l Jv o,ip + 9*1 
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This is the case for the A = 0.06 and A = 0.07 models shown in the left panel of Fig. [4] These bounds are independent 
of the decay rate T x . Furthermore, the bound Eq. (37) is written solely in terms quantities evaluated at horizon 
crossing, and hence may be computed without explicitly knowing the full non-linear evolution of the bundle during 
the reheating process. Whether this maximum value, \f^ Imaxj is obtained at the end of reheating is of course 
dependent on T x . Formally, the lower bound for /j\jL a ' (which is approached as T x — > 0) would be /NL, al m in = 0- 

The existence of a minimum of Eq. (35) at negative N\ tV> for A = 0.05 explains the seemingly opposite dependence 
of /^L al on T x compared to A = 0.06 where the minimum exists at positive Niu>: as T x is increased from zero (where 
Ni -ip — > — oo), the time taken for reheating to complete is reduced and tp freezes out sooner, hence reducing the 
mag nitu de of N\ tV . For the A = 0.05 model, as T x is increased further, driving N\, v toward zero, the minimum of 
Eq. (35) is encountered, past which point I/nl^I ^ s reduced. For A = 0.06, increasing T x still drives N% >{ p toward zero, 
but this time |/^ al | is increased. 

For A = 0.07, the function Eq. (35) is almost completely flat for Ni >(p < 0, which indicates that no matter how 
slowly or rapidly the universe is reheated, the value of /nl at the start of reheating will survive until completion. In 
the limit of instantaneous reheating, T x — > oo, Ni -ip sa 0, and so f^f^ ~ /nl(M- This is only approximate since, as 
reheating does not begin on a hypersurface of constant density, there will be some small correction N% v . 

Another interesting observation is that I/nl^I ( or m ore accurately the derivative Ni tV ) is fairly insensitive to 
changing the decay rate by many orders of magnitude. For example, as can be seen from Table [ij I/nl 3,1 ! changes by 
less than three units as the decay rate is increased from T x = \/lO~ 5 WoM p to T x = y/10~~ 1 WoM p . We caution here 
that decay rate could, in principle, be many orders of magnitude weaker than the weakest decay rate studied here 
and still be consistent with the bound derived from BBN constraints, T x > 4 x 10 _40 M p . These tiny (but non-zero) 
values of T x are beyond our numerical capabilities: to compute the statistics of £ at the completion of reheating 
requires integrating the field equations up until the universe is radiation dominated, which for such weak rates, can 
take O(30) e-folds. Substantial errors are accumulated if the field equations are integrated over such long periods of 
time, which in turn induces large errors in the computation of the SN derivatives. For this reason, we only quote 
values of /nl, an d r f° r decay rates for which we are confident that we have control over all sources of numerical 
error. 



B. Quartic minimum: W(tp, %) = WoX 4 e 



We now repeat the same analysis, but with a quartic minimum in the x direction. The background inflationary 
dynamics are similar to the X 2e ~ model as can be seen from the slow-roll solutions to the field equations: 

X 2 = xl-8N, ^ = ^e 2AAr . (39) 

The oscillatory dynamics about the minimum are somewhat different to that of the x 2 case however, due to the 
potential being much shallower around x = 0, with steep sides away from the minimum. The oscillations of x are not 
sinusoidal, but are instead given approximately by the elliptic function [58, 66 when T x = 0: 

* T >"l/s*i3F "(?■•£)• <40) 

where r is the conformal time dt = a(t)dr. Here, b ss 0.85 is a numerical constant and u> is the effective frequency 
of the oscillations. The energy density of the coherently oscillating x field decreases in the same way as a relativistic 
fluid, p x <~~> a~ 4 , behaving as radiation with a non-vanishing pressure P x « \Px- 

Provided A is not too large, the tp field remains slowly rolling throughout the entire reheating phase. In the left 
and right panels of Fig. [5] we show the final stages in the evolution of /nl and N v respectively as a function of N 
for various decay rates T x . We see that the qualitative dependence of /n'l*' 011 the decay rate is the same as for 



the x e model, which may be explained by appealing to Eq. (35). This implies that the shape of the minimum 



does not change the qualitative dependence of /nl 11 011 the reheating process. Of course, as reheating proceeds, the 
shape of the x minimum does not remain exactly quartic (or quadratic in the case of the previous model) due to the 



coupling with the tp field. We will comment more on this in Section VI 



C. Spectral index and tensor— to— scalar ratio 



For completeness, we also examine how sensitive the tensor-to-scalar ratio r, and spectral index n^, are to the 
reheating phase for the models studied above. From Eqs. (18) and (20), we have that r = S/M p Q2i N /)• Recall 
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FIG. 5: W(<p,x) = Wox 4 e~ Xv . The parameters used are: A = 0.055, ip, = 5 x 10" 4 M P and x * = 23.0M P . Left Panel: The 
evolution of /nl during reheating for various decay rates F x . Right Panel: The evolution of the derivative N v during reheating 
for various decay rates F x . All decay rates are in units of y/WoM p . In both panels, the solid vertical (black) line denotes the 
end of inflation, 7V C , and the dashed vertical (blue) line denotes the start of reheating, iV r . The Hubble rate at the start of 
reheating is H r « ViO -1 WoM p . 



that due to the hierarchy in magnitude between the scalar field kinetic energies ipq and \* a t horizon exit, we can 
approximate <?* = N x w (2e* ) -1 / 2 . For the region of parameter space of interest, as T x is decreased from infinity, the 
time taken for reheating to complete is increased and ip freezes out later, increasing the magnitude of Ni lV . Hence, 
the weaker the decay rate, the more suppressed the tensor-to-scalar ratio, and the following bound exists: 



r < 



1 



M 2 N* 



Q,<p + .9* 



(41) 



This suppression of r for weaker T x is illustrated in Table [T] 

A similar bound also exists for the spectral index. Whilst it is a good approximation to neglect N vx in the expression 
for /nl> one must be more careful in the expression for uq: the hierarchy -C \x*\ means that the sum in the 
numerator of Eq. (21 1 generates terms of similar order. 



X 2 minimum: /nl(^) = —5.93, 
n ( {t e ) = 0.763, r(i e ) = 2.8 x lfT 4 



X 4 minimum: fm,(t e ) = -48.29, 
n s (t e ) = 0.770, r(t e ) = 7.2 x 10~ 3 





£ final 

/nl 


^final 
"g 


™final 




r x 


f final 

/nl 


„ final 
"s 


„ final 


\/l0- 5 


-4.35 


0.761 


2.4 x 10" 


-4 


Vio- 8 


-54.40 


0.772 


9.7 x 10" 3 


V10- 3 


-5.54 


0.762 


3.9 x 10" 


-4 


Vio- 6 


-60.32 


0.778 


1.2 x 10" 2 


Vio- 1 


-7.14 


0.762 


6.3 x 10" 


-4 


Vio- 4 


-65.80 


0.776 


1.5 x 10" 2 



TABLE I: Statistics of £ for W(<p, X ) = WoX a e Xlf2 for different decay rates. All decay rates are in units of y/WoM p . We 



give values computed at the end of inflation (t e ) 
Quadratic minimum (a = 2); A = 0.06, ip* = 10" 
tp* = 5 x 10~ 4 M P and X * = 23.0M P . 



and at the completion of reheating (final) where £ is conserved. Left Table: 
3 M P and x* = 16.0M P . Right Table: Quartic minimum (a = 4); A = 0.055, 



However, the expression for nt, Eq. (21), can be reduced to a function of solely N v , by making use of the scaling 
relations Eqs. (p9l and ffl. Using g* = N x w (26*)" 1 / 2 and N xx w 1 - (r) xx /2e x )« = 1/2 in Eq. (Ell, with 
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(<Pi/H). 



-(2e* I ) 1 ^ 2 and relations Eqs. (29) and (30) we may write: 



-2e, 



N 2 

•p 



4A 




(42) 



By the completion of reheating, the last two terms in parenthesis are small, suppressed by factors of N v and may be 
neglected. With J2e* ftp* — 2A we arrive at 



nc - 1 -2e* 



_4XNl_ 
N 2 +g; 



> -2e* -4A. 



(43) 



This shows that there is a bound for the allowed range of uq, and also explains the almost complete insensitivity of 
tiq to T x when N 2 3> g 2 , an example of which is given in Table |l] for A = 0.06. By comparison, /nl is much more 
sensitive to T x for the same value of A, since /nl(^) is n °t na t over the N v range of interest. This indicates that, 
for this particular inflationary model, tiq is a more robust inflationary observable, and perhaps a better probe of 
the underlying inflationary model since it is insensitive to the physics of reheating. Whilst the two models (with 
quadratic or quartic minima) studied in this section can generate a large /nl that survives until the completion of 
reheating, they would be ruled out by observation since their spectral indices are far too low. 

The qualitative arguments given in this section are respected as long as A is not too large. If ip ~ 0(A -1 / 2 ), then the 
scaling relations Eqs. (29) and (30) will break down and our arguments may no longer hold. The effect of reheating 
on the motion of ip is to prevent it from rolling any further down its potential: as reheating becomes more efficient, 
ip freezes out at smaller values. In this strong coupling regime, we expect the scaling relations to work well, but 
becoming a worse approximation if reheating proceeds very slowly. 



IV. TWO MINIMA 



Assisted inflation [99, may be realised via a collection of string axions. In this scenario, known as N-fiation [100J, 
the many axion fields cooperatively source inflation even if their potentials are individually too steep. The collective 
potential is comprised of a sum of Nf uncoupled axions <p>c 



N f 

W(cp) = '£A* 



1 — cos 



2tt 

li 



■Pi 



(44) 



With only a single field present, this model is more commonly known as natural inflation |101j . Each axion is fully 
described by its decay constant fi and its potential energy scale A 4 . The standard arguments show that we should 
expect fi ~ 10 16 GeV. The mass of each field in vacuum satisfies n^t,u\ — 47r 2 A 4 // 4 2 . Due to the shift symmetry 
<Pi n fit we can without loss of generality set the initial conditions tp*u) <E [0, fi]. To generate a large /nl, we 

must have at least one axion close to the 'hilltop' at p ~ 0.5 [102) . The evolution of the curvature perturbation in the 
post-inflationary epoch, including the reheating stage, for sum-separable multifield models was also studied by Choi 
et.al. 



[103J. However, these authors restricted themselves to scalar fields with quadratic minima and similar masses 
which decay with similar rates; while in the following we consider a more general scenario and include the studies of 
non-gaussianity as well as the tensor-to-scalar ratio r. 



A. Quadratic minimum 



We follow [38] , supposing that the initial conditions are chosen so that only a single axion ip populates this hilltop 
region. This field sources the non-Gaussianity, whilst the remaining Nf — 1 axions, which begin away from the hilltop, 
contribute only to the expansion rate. By expanding about the minimum of the remaining Nf — 1 fields, these axions 
may be replaced by a single effective field \ with a quadratic potential. With fi = f for all axions, the effective 
two-field potential then reads: 



W(<p,x) = W 



^mV + A 4 ^1 - cos 



(45) 
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FIG. 6: W{if,x) = W [|m 2 x 2 + A 4 (l - cos (yv))]- The parameters used are: A 4 = rrff/A-n 2 , <p, = (| - 0.001)M P , 
X* = 16A/ P , f — m — 1. Both panels show the evolution of /nl during reheating. Left Panel: Equal decay rates, T x =Y ^ ^ 0. 
For comparison we also show the Y x = Y v = limit (thin black line). Right Panel: Unequal decay rates, r x 7^ Y v 7^ 0. For 
comparison we also show the Y x = Y v — limit (thin black line). In both panels, the solid vertical (black) line denotes the 
end of inflation, N c , the dashed vertical (blue) line denotes the start of x reheating and the dotted vertical (red) line denotes 
the start of <p reheating. The background Hubble rates at the x an d f reheating surfaces are H* « y/5 x 10~ 2 WoM p and 
Hf w VW- 2 W M P respectively. 



In fact, replacing the collective potential with an effective two-field potential is well motivated, see for example [104J, 
where they showed that the energy density of the universe is dominated by fields with comparable masses even if 
one starts with thousands of fields, including the post-inflationary reheating stage. Reheating in models of N-flation 
also proceeds preferentially via a perturbative decay route as opposed to via parametric resonance and preheating 

From this point onwards we will refer to ip as the axion and to x as the inflaton. By suitably choosing the 
axion/inflaton mass ratio in vacuum, various scenarios can be realised. For example, if the axion is sufficiently 
massive it may quickly decay to its minimum during inflation, where it becomes trapped without oscillating. In this 
case, adiabaticity is established long before reheating begins, and the decay of the inflaton into radiation does not 
affect the evolution of £. We have confirmed this numerically. 

It is also possible to realise dynamics where both fields minimise after inflation has ended, entering an oscillating 
phase such that perturbative reheating can be applied. For example, with A 4 = m 2 / 2 /47r 2 , (p* = (| - 0.001)M p , 
X* = 16Mp and / = m = 1, the inflaton minimises before the axion, but both fields minimise after inflation has 
ended. In this example both fields acquire the same mass in vacuum. Fig. [6] shows the evolution of /nl for different 
combinations of T x and T v for this parameter choice. Since both fields oscillate rapidly about their minima, both 
fields must be coupled to radiation. If one field is instead left uncoupled, its energy density will scale as matter since 
the minimum is quadratic, and will eventually come to dominate over radiation which redshifts away more quickly. 

Unlike the product separable case where the universe is reheated from only a single field, the ip field has left slow-roll 
by the time reheating starts. Hence, the non-linear dynamics during the oscillating phase is essential and we could 
not find any simple scaling relation between N vtp , N vx and N v . Yet we find that /nl is still dominated by the same 
term as in the case where adiabaticity is reached before inflation ends [38] : 



5JV, 

/NL ~ 



(46) 



6 TV 2 ■ 

As can be seen from the left panel of Fig. [6j /nl 3,1 i s almost completely insensitive to reheating when T x ~ T v . 
However, as can be seen from the right panel, a mild hierarchy between T x and T v generates significant corrections to 
to /nl, 3 '- This effect is not due to the axion reheating hypersurface being distinctly separated from the inflaton surface 
(the vertical dotted (red) and dashed (blue) lines of Fig. [^respectively) and we have confirmed this numerically. What 
is important however, is the axion/inflation mass ratio in vacuum. The model parameters which realise the dynamics 
seen in Fig. [6] give m v = m x at the minimum. The differences induced in /^L al when a mild hierarchy exists between 



17 



r x and is greatest when the masses are equal. As the masses are separated, keeping the ratio T x /T ip fixed, the 
sensitivity of f^1 al to reheating decreases. This can be understood as follows: first consider the situation where the 
two fields have different masses, for instance, m x > m v . Assuming both fields reheat at roughly the same time, 
the more massive field \ wm dominate the energy density of the universe and thus the dynamics of the universe 
during reheating. Evaluating on constant energy hypersurfaces, the initial horizon crossing dependence of the x field 
dynamics is smaller compared to the case m x = m v , where the energy density of the universe is distributed evenly 
between the fields. As a result, we expect the number of e-folds of expansion N and are less sensitive in the 

case m x ^ m v . 

In fact, having the two fields decay at different rates is a form of modulated reheating, although it is different from 
the standard scenario 18-20]. In the standard modulated reheating scenario, inflation is driven by a single field, 
whose decay rate is modulated by a second, subdominant field that remains light and plays a negligible role during 
inflation. The fluctuations of the subdominate field induce fluctuations in the inflaton decay rate and thus generate 
curvature perturbation during reheating. In the two minima case here, note that the initial horizon crossing values 
of the fields tp*,x* determine how the energy density of the universe is distributed between the two scalar fields. 
Therefore, although the field decay rates are constant here, the rate of energy transfer from the scalar fields to the 
radiation fluid can be different for each inflationary trajectory in the bundle and thus can generate extra contributions 
to the curvature perturbation, provided there is a mild hierarchy in the decay rates. Therefore it is not surprising 
that /nl can acquire such a significant correction during reheating when the two decay rates are different. The 
two minima scenario is also similar in spirit to a model of two field inflation with equal masses followed by instant 
preheating, in which the two field have very different couplings to the preheat field [106] . for a related scenario see 
also [107] . Note however that all of these instant preheating models are very tightly constrained even at the level of 
linear perturbations |21j . 



B. Quartic minimum 



We now repeat the same analysis, promoting the quadratic \ 2 minimum of Eq. (45 1 to a quartic minimum, x 4 . This 
modification was also studied in [108] where the model parameters were chosen such that £ becomes conserved during 
slow-roll. Again, we find a similar qualitative behaviour of /Ifff 1 as in the quadratic case: the asymptotic values 
of /nl are very insensitive to the decay rates of the scalar fields when they are equal, and slightly more sensitive if 
they are different. However, all observables are much less sensitive to decay rates here as compared to the quadratic 
minimum case. 

We summarize the values of the observables of £ for the models studied in Section IV A and IV B at the end of 
inflation and end of reheating in Table [TTl 



X 2 minimum: /Nx(*e) ~ 0, % 4 minimum: /nl(^) ~ 0, 

n s (t e ) = 0.969, r(Q = 0.124 n s (t e ) = 0.951, r(t e ) = 0.263 
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TABLE II: Statistics of C for W(<p, X ) = W [\m 2 x a + A 4 (l - cos (^j J J for different decay rates. All decay rates are in 

units of y/WoMp. We give values computed at the end of inflation (t e ) and at the completion of reheating (final) where C, is 
conserved. Left Table: Quadratic minimum (a = 2); A 4 = m 2 f 2 /4iT 2 , ip* = (i — 0.001)Af p , \* ~ 16M P , / = m = 1. Right 
Table: Quartic minimum (a = 4); A 4 = m 2 f / 4tt 2 , tp* = (| - 0.001)M P , \* = 22M P , / = m = 1. Notice the very large 
decrease in the tensor-to-scalar ratio from the end of inflation to its final value. 
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FIG. 7: W(x,<P) = W (x 4 e- Xv /M p + k X 2 )- We show / NL as a function of N during reheating. The parameters used are: 
ip* — 10 _3 M P , x* = 22M P and A = 0.06. In both panels, the solid vertical (black) line denotes the end of inflation, N e , and the 
dashed vertical (blue) line denotes the start of reheating, N T . Left Panel: k = 1.0. The Hubble rate at the start of reheating 
is H r « V2 x 10 "WoMp. Right Panel: k = 0.1. The Hubble rate at the start of reheating is H r w VlO -1 W M P . 



V. NON-SEPARABLE POTENTIAL WITH ONE MINIMUM 



In previous sections, we have studied the evolution of /nl and its asymptotic value at the end of reheating, 
in examples where one or both fields reheat from a two-field separable potential. In this section, we will repeat the 
same analysis, but this time for a non-separable potential. 

As an example, we consider a modified version of the previously studied quartic exponential model, by adding an 
extra quadratic mass term 

W( X , <p) = W (x 4 e-W M p + K X 2 ) . (47) 

Before discussing reheating, it is useful to first study the inflationary regime. During inflation, the quadratic x 2 mass 
term has a negligible effect on the field dynamics when the x field is of O(l) in Planckian units, unless k 3> 0(1) or 
Xip 2 3> O(Mp). Here in the following, we only consider the case k ~ 0(1), for which we can approximate the field 
dynamics and /nl during inflation as the same as setting k = 0. Therefore, in the region of parameter space where 



k < O(l), /nl is expected to follow similar evolution as in the separable case studied in Section IIIB during the 
slow-roll regime, with large deviations only coming in at late times towards the end of inflation. 

The mechanism for generating large /nl is the same as discussed in [35] , which is well illustrated from the fact that 
there exists a scaling relation for the subdominate field 5N derivatives. 

For the values k = 1, A = 0.05, y>* = 10~ 3 Af p , x* — 22M P , a large negative /nl is generated during inflation as 
the (p field rolls down the ridge and the bundle of trajectories diverge. The evolution is similar to the separable case 
where K = 0, with /nl ~ —44 close to the end of slow-roll. Things are however a bit different after inflation even 
before reheating starts. For A = 0.06, the additional quadratic term becomes comparable to the quartic term slightly 
earlier than in the case of A = 0.05. In this case, we find /nl swaps sign shortly after the end of inflation. This 
unexpected behaviour, which we do not see in other cases, may be explained as follows: although the trajectories 
are still diverging in the <p direction in this case, the fact that the quadratic term becomes dominant suggests that 
the local potential geometries around each trajectory converge to the same quadratic shape, independent of ip. This 
would have the same effect as the trajectories themselves converging in the separable case where H is converging, 
thus giving momentarily large positive /nl- 

Shortly after inflation ends, when the x ne ld reaches sub-Planckian values, the x 2 term starts to dominate over 
the x 4 term. Therefore, we expect the additional % 2 term modifies the field dynamics during the reheating phase and 
possibly /nl as well. The additional x 2 term makes the potential less shallow around the minimum. This saves the 
X field from being frozen to non-zero values, leaving unwanted residual potential energy in the case where T x is too 
large where the oscillations of the scalar fields are heavily damped. This happens if the potential around the minimum 
is too shallow, as in the model studied in Section |IIIB| Similar to the separable case, as shown in Fig. [7J we found 
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during the early reheating stage, with a larger amplitude for smaller T x . 

are now much less sensitive 



that /nl oscillates roughly in phase with x 2 
However, unlike the previous separable case in Section III B 
to T x and thus the reheating timescale. The sensitivity increases as k decreases as shown in Table III This means 
the effect of introducing additional quadratic mass term reduces the sensitivity of /nl to the reheating timescale. 6 



the 5N derivatives and /] 



NL 



Non-separable n = 1.0 /nl(^c) = —18.71, Non-separable k = 0.1 /nl(^c) = —13.23, 

n s (t e ) = 0.748, r(i e ) = 4.1 x 10~ 3 n s (t e ) = 0.746, r(t e ) = 2.0 x 10~ 3 
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TABLE III: Statistics of ( for W(x, f) = W (x i e~ X ' p /Mp + ^x 2 ) for different decay rates. All decay rates are in units of 
V Wo Alp. We give values computed at the end of inflation (t e ) and at the completion of reheating (final) where ( is conserved. 
Left Table: The parameters used are: (p* = 10 _3 M P , \* = 22M P and A = 0.06 and k = 1.0. Right Table: The parameters used 
are: ip* = 10" 3 M P , x* = 22Af p and A = 0.06 and k = 0.1. 



VI. DISCUSSION AND CONCLUSIONS 



In this work we have studied the effects of perturbative reheating on the key inflationary observables /nl, an< ^ r ' 
for canonical two-field inflation models. We have considered two classes of potential: the 'runaway' type which has a 
minimum in only one direction; and potentials which have a minimum in both directions. We have studied quadratic 
and quartic minima, finding that the dependence of the statistics of ( on the decay rate of the field (s) is qualitatively 
the same. Perhaps most importantly, we have shown for both classes of models, that if a large non-Gaussian signal 
exists at the start of reheating, it will in general be non-zero at the completion of reheating. 

For the single minimum models, adiabaticity is never established before reheating begins and so the bi-spectrum 
acquires substantial reheating-dependent corrections. As a consequence, the magnitude of any non-Gaussianity 
generated at the end of inflation does not necessarily remain the same at the end of reheating, meaning that /nl 
cannot be linked directly to the physics of the inflationary model. Whilst /nl is sensitive to reheating, we have also 
shown that there can exist certain regimes of model parameter space where the spectral index is almost completely 
insensitive to reheating. In such scenarios, may be considered a more robust inflationary statistic and a better 
probe of the underlying potential. 

For two-minima models where both fields decay to reheat the universe, we have shown numerically that even if 
an adiabatic condition is approached by the inflating/isocurvature fields converging in, and oscillating about, their 
global minima, the decay of these fields into radiation can promote further evolution of £. If a mild hierarchy in decay 
rates between each field exists, /nl can be enhanced or suppressed relative to the same model where reheating is not 
accounted for. 

One important difference between the single minimum models and the two minima model is that in the former 
case, the fields are coupled via the potential, whilst in the latter they are coupled only via gravity. Thus, for the 
single minimum models of Sections III and |Vj the local geometries of the x minima are functions of the sub-dominate 
field tp, and these geometries are different for different inflationary trajectories in the bundle. This is illustrated in 
Fig. [8j The shape of these 'reheating minima' evolve in time as reheating proceeds, and will affect the dynamics of 
the oscillating x field. In the two-minima model of Section |TV| however, where the potential is sum-separable and the 
fields are coupled only through gravity, the local geometries of the x minima are always independent of tp and so this 
effect is not present. This effect is of course model dependent, as we have illustrated with the non-separable model 
Eq. (47). When the interaction term is small and plays a negligible role during reheating, the x held dynamics are 
independent of the dynamics of ip. This explains why we found the sensitivity of the SN derivatives to T x decreases 



Note that changing k also slightly changes the times that inflation ends and reheating starts. This would however have negligible effect 
on the dependence of the observables on T x in the parameter space of interest. 
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FIG. 8: An illustration of how the reheating minima depend on the subdominate field as in the case of quadratic times 
exponential potential. The perturbation in the <f> direction is denoted by S. Left Panel: Local geometries of the minima around 
the x direction. Right Panel: The potential gradient in the x direction. 

as k increases. For larger k, the ip dependence of the local geometries of the x minima decreases and thus /nl is 
very insensitive to T x . 

In general, for both classes of potential (one and two minimum), the degree of sensitivity of the primordial observ- 
ables to reheating seems to be a model (and parameter) dependent issue. That said, for all the models studied in this 
work, whilst the magnitude of /nl may be heavily dependent on the decay rate of the field(s), its sign remains the 
same. Whether this conclusion encompasses more complicated models is unclear. For example, couplings between 
the inflating and isocurvature fields may promote the global minima studied in this paper, to functions of the fields 
themselves (local minima) possibly leading to a highly non-trivial reheating surface which may generate large cor- 
rections to £. Models with non-canonical kinetic terms, such as DBI inflation may also not respect this observation, 
since the fields no longer follow the gradient of the scalar potential. 

Another point worthy of some discussion is that we have assumed reheating to take place entirely perturbatively. 
Whilst we consider this to be a sensible starting point for a first-time exploration of the sensitivity of /nl, ano - 
r to reheating, this is almost certainly a gross simplification: the initial stages of particle production (preheating) is 
a violently explosive non-perturbative effect. It is expected that such a rapid preheat stage in the regime of broad 
resonance may have long-lasting effects on the subsequent evolution of the universe. For example, it may lead to 
specific non-thermal phase transitions in the early universe, (1091 II l(Jj . topological defect production and promote 
novel mechanisms for baryogenesis [1111 1112] . Preheating has also been shown to generate significant levels of non- 
Gaussianity O HH HM2Tj . If a more sophisticated description of (p)reheating were employed, including the rich 
spectrum of perturbative and non-perturbative QFT effects, it is tempting to speculate that the statistics of £ might 
be more radically altered in models of multi-field inflation. 

Non-Gaussianity has evolved into a very active and topical field, in which observations have improved greatly over 
the last decade, through both studies of the CMB and large scale structure. At the same time, on a theoretical 
and phenomenological level, a plethora of different mechanisms have been suggested which are capable of generating 
an observable /nl- Currently, the tightest constraints on local type /nl come from the WMAP satellite, which 
constrains the amplitude of the non-Gaussian part of £ to be less than about one thousandth of the amplitude of 
the Gaussian perturbation. Planck is expected to tighten this constraint considerably. A detection of /nl at this 
level would rule out the simplest canonical, single field inflation models, where it has recently been explicitly shown 
that preheating has a negligible effect on the scalar bi-spectrum [52]. However, we will most likely be left with many 
other viable scenarios, including multi-field models, which when suitably tuned, can match the observations. As 
we have demonstrated in this paper, accounting for the dynamics of reheating muddies the waters further. Unless 
adiabaticity has been achieved before the onset of reheating, it seems unlikely that we can use explicit values of /nl to 
discriminate between different multi-field models unless we have a complete understanding of the reheating process. 7 



This applies to local— types of non-gaussianity only. For other shapes such as equilateral type, the post-inflationary evolution will not 
change the inflationary predictions as the contributions come from interactions at or before horizon-crossing, see for examples ITT 3 115 
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In this sense, our work lends support to that of |116j . in that it also represents a challenge to the conventional lore 
that non-Gaussianity is a 'smoking gun' signature of non-standard inflationary dynamics: such signatures may be 
significantly altered by the subsequent reheating phase. 

More optimistically, non-Gaussianity is not only about one single number. The trispectrum (the four-point func- 
tion) depends on two non-linearity parameters tnl and <?nl, and if the current observations |/nl| ~ 40 (which are 
not statistically significant) turn out to be true, then tnl should be large enough for Planck to detect. It might be 
that, like the spectral index for the single minimum models studied in this paper, the trispectrum is less sensitive 
to the physics of reheating. Furthermore, if /nl is detected, it may also be possible to constrain or even detect a 
scale dependence: /nl is often assumed to be constant, but this is only true for certain simple models. For example 
/nl is strongly scale dependent in the two-field hybrid inflation model [24]. This opens up the question of whether 
reheating may leave some observable signature in the running of /nlj which may be used as a complimentary probe 
of the inflationary theory and the reheating mechanism itself. Indeed, whilst it has been shown that /nl is insen- 
sitive to preheating in canonical single field models (as well as being too small to be observed) it is strongly scale 
dependent [52] , 

When discussing the sensitivity of the primordial observables to reheating, it is also important to keep in mind the 
degree of fine tuning which is required for the inflationary model itself to be consistent with observations. For typical 
models this amount of fine tuning is large, especially if one wishes to generate an observable /nl- Accounting for 
the subsequent dynamics of reheating introduces a further source of fine tuning, however what is apparent from this 
work, is that this is secondary compared to the degree of inflationary fine tuning. For the single minimum models 
studied in this paper for example, changing A or ip* by one part in 10 2 may completely remove any non-Gaussian 
signal, whilst shifting the reheating decay rate by two orders of magnitude changes /nl by 0(2) units. 

In conclusion, whilst non-Gaussianity is in principle a powerful probe that may be used to distinguish between the 
many models of inflation, we must be careful in our interpretation of any observational constraints that place bounds 
on the statistics of £. We have shown that the dynamics of perturbative reheating can have a non-negligible impact 
on these statistics for canonical two-field inflation models. As such, without a UV complete theory of inflation and 
reheating, it seems hard to infer the properties of the underlying inflationary potential from observational bounds on 
/nl and related quantities alone. 
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